      !<@file

      subroutine inversao()

      use globais
      implicit none


      character*5 :: filename1,filename2,filename3
      double precision, dimension(:), allocatable :: X !< variaveis de otimizacao/parametros da inversao
      double precision, dimension(:), allocatable :: X0 !< variaveis de otimizacao/parametros da inversao iteração anterior
      double precision:: cmin0,cmax0
      double precision:: xdonoh,zdonoh
      integer :: i,j
      double precision :: VALUE	
      integer :: kncl
      integer :: Nclasses !< Melhor número de classes para clusterização.
      double precision :: mapCemX
      double precision:: tolerro=0.02d0 !< maximo erro admissivel para a razão entre a norma do residuo e a norma do vetor medidas

      real, dimension(:,:), allocatable:: creal !< Matriz auxiliar, velocidade de propagação em real*4 para impressão em arquivo binário a ser plotado pelas rotinas do Seismic Unix.

      integer :: MaxIter
      integer :: kregaux

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
! ! non linear conjugate gradient, napack, netlib
      !input
      double precision :: STEP  !--STEP SIZE ALONG SEARCH DIRECTION FOR FINAL ITERATION		   
      double precision :: Tol  !--COMPUTING TOLERANCE (ITERATIONS STOP WHEN MAX-NORM OF GRADIENT .LE. Tol)  
      integer :: limit
      double precision , dimension(:,:), allocatable::Work
      !output
      double precision :: E  !--MAX-NORM OF GRADIENT


      EXTERNAL BOTH,GRAD,PRE,VALUE
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 


      !< L-BFGS 
!       PARAMETER(NDIM=2000,MSAVE=7,NWORK=NDIM*(2*MSAVE +1)+2*MSAVE)
      DOUBLE PRECISION, dimension(:), allocatable:: DIAG
      DOUBLE PRECISION, dimension(:), allocatable:: W
!       DOUBLE PRECISION, dimension(:), allocatable:: G
!       DOUBLE PRECISION, dimension(:), allocatable:: G1
!       DOUBLE PRECISION, dimension(:), allocatable:: G2

      DOUBLE PRECISION, dimension(:), allocatable:: Qlbfgs
      DOUBLE PRECISION, dimension(:), allocatable:: Q1lbfgs
      DOUBLE PRECISION, dimension(:), allocatable:: Q2lbfgs
      DOUBLE PRECISION, dimension(:), allocatable:: RHOlbfgs
      DOUBLE PRECISION, dimension(:), allocatable:: ALPHAlbfgs
      DOUBLE PRECISION, dimension(:), allocatable:: Slbfgs
      DOUBLE PRECISION, dimension(:), allocatable:: Ylbfgs
      DOUBLE PRECISION, dimension(:), allocatable:: Y1lbfgs
      DOUBLE PRECISION, dimension(:), allocatable:: Y2lbfgs


      DOUBLE PRECISION :: EPS,XTOL,GTOL,STPMIN,STPMAX
      INTEGER ::  IPRINT(2),IFLAGlbfgs,ICALL,MP,LP
      LOGICAL ::  DIAGCO
      integer :: msave
      integer :: icallmax
      integer :: itertotal
      DOUBLE PRECISION :: FobjetivOLD!<valor da função objetivo na avaliação anterior, usado como referência para critério de parada
      DOUBLE PRECISION :: MisfitOLD!<valor da função residuo na avaliação anterior, usado como referência para critério de parada
      DOUBLE PRECISION, dimension(:), allocatable:: bestX !< melhor X até o momento
      DOUBLE PRECISION :: bestFobjetivo !< valor da função objetivo correspondente ao bestX      
      DOUBLE PRECISION :: bestMisfit !< valor do resíduo correspondente ao bestX      

  


      EXTERNAL LB2
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 

      double precision:: fobjetivo !< funcao objetivo


      write(*,*)'Solver Inversao'


      call entradainversao(maxIter,DX,tol,limit,eps,xtol,msave,cmin0,cmax0)


      write(*,*)'alocando X0...'
      allocate(x0(n))

      write(*,*)'alocando DX...'
      allocate(DX(n))


      write(*,*)'alocando X...'
      allocate(x(n))


      write(*,*)'alocando bestX...'
      allocate(bestX(n))


      write(*,*)'alocando cX...'
      allocate(CX(n))



      

      do j=1,nz0
	  do i=1,nx0
	      X((j-1)*nx0+i) = cmin0  + (cmax0-cmin0)/dfloat(nz0-1)*(j-1) !< aproximação inicial linear da velocidade. Falta transformar no parâmetro X
	  enddo
      enddo


      do i=1,N
	  X(i) = mapCemX(X(i),cmin,cmax) !< aqui, X(i) argumento ainda é velocidade. Volta da função o parâmetro X correspondente
      enddo




      filename3= 'sol_c'
      call system('rm '//filename3//'.bin '//filename3//'.ps')
      open(unitsol,file=filename3//'.bin',status="unknown",form="unformatted",access="direct",recl=4*Nx0*Nz0)      

      open(unitF,file='outfuncao.txt',status="unknown")
      write(unitF,*)'1: freq em Hz, 2: iter, 3: misfit, 4: lambda, 5: treg, 6: misfit relativo (%)'
      open(unitFgrad,file='outfuncao_grad.txt',status="unknown") 
      write(unitFgrad,*)'1: freq em Hz, 2: iter, 3: misfit, 4: lambda, 5: treg, 6: misfit relativo (%)'

      write(*,*)'alocando Xn...'
      allocate(Xn(n))
      write(*,*)'alocando cxn...'
      allocate(CXn(n))
      Xn = X


      if(metodo.eq.2)then
	  open(unitBFGSlp,file='BFGSlp.txt')
      	  open(unitBFGSmp,file='BFGSmp.txt')
      else
	  call system('rm BFGS*.txt')
      endif

      kregaux=kreg! salva o carreg original

      open(unitX,file='outX.txt',status="unknown")

      do kloop=1,Nloops
  
! 	  if(kloop.gt.1)kreg=0 ! sem regularização apś a 1a rodada nas frequências 
! 	  kreg = kregaux
! 	  if(mod(kloop,2).eq.0)then ! sem regularização nas rodadas pares
! 	      kreg=0
! 	  endif
	  !do kw=kloop,Nomega,Nloops          !<loop seletivo, seleciona freqs de acordo com o numero de loops
	  do kw=1,Nomega

	      
	      write(unitF,*)! só no inicio das novas frequencias,salta linha
	      write(unitFgrad,*)

	      MisfitOLD = 0.0d0
	      fobjetivOLD = 0.0d0
	      itertotal=0
	      lambdanaux=1.0d0
	      bestFobjetivo = 10d20
		
! 	      if(kw.eq.2 .or. kw.eq.11 .or. kw.eq.19 )then!seleção das freqs
	      write(*,*) '#menor comprimento de onda (pelo cmin)= ', cmin/(omegas(kw)/(2*Pi))
	      write(*,*) '#intervalos por comprimento de onda (pelo cmin)= ',  ( cmin/(omegas(kw)/(2*Pi)))/delta
	      write(*,*)' kw ',kw,': ',omegas(kw)/(2.d0*Pi),'Hz'

	      norm2PBarra=0.d0
	      do igsr=1, gruposSR
		  do ks=1,Nsource(igsr)
		      do kr=1,Nrec(igsr)
			  norm2PBarra= norm2PBarra + PreceivBARRA(kr,ks,kw,igsr)*CONJG(PreceivBARRA(kr,ks,kw,igsr))
		      enddo
		  enddo
	      enddo
	      write(*,*)'||Pbarra||**2=',norm2PBarra

	      if(metodo.eq.1)then

		  write(*,*)'alocando work para o cg do napack...'
		  allocate(Work(N,3))	
		  iterant = -1! numero da iteracao anterior
		  write(*,*) 'step:',step,'tol:',tol,'iteracoes por normalizacao', maxiter,'max. iter do CG:',limit
		  call CG(X,E,ITER,STEP,Tol,LIMIT,N,MaxIter,VALUE,GRAD,BOTH,PRE,Work)
		  write(*,*)'fim do CG   step:',step
	! 		    step=0.0d0
		  deallocate(Work)! do cg do napack

	      elseif(metodo.eq.2)then

		  allocate(DIAG(N))
		  allocate(W(N*(2*MSAVE +1)+2*MSAVE))
! 		  allocate(G(N))
! 		  allocate(G1(N))
! 		  allocate(G2(N))
		  allocate(Qlbfgs(N))
		  allocate(Q1lbfgs(N))
		  allocate(Q2lbfgs(N))
		  allocate(RHOlbfgs(msave))
		  allocate(ALPHAlbfgs(msave))
		  allocate(Slbfgs(N*msave))
		  allocate(Ylbfgs(N*msave))
		  allocate(Y1lbfgs(N*msave))
		  allocate(Y2lbfgs(N*msave))

 15		  ICALL=0
		  maxiterLBFGS = limit
		  icallmax= 2*maxiterLBFGS
		  IPRINT(1)= 1
		  IPRINT(2)= 0
		  DIAGCO= .FALSE.

		  iterant = -1! numero da iteracao anterior
		  iter =0


		  do i=1,n
		      X0(i) = 0.0d0
		  enddo


		  IFLAGlbfgs=0
		  
		  write(*,*) 'eps:',eps,'xtol:',xtol,'maxiterLBFGS:',maxiterLBFGS
		  DO
		      write(*,*)
		      write(*,*)
		      write(*,*)
		      write(*,*)'---------------------------- LBFGS:',iter,'(iter)',itertotal,'(itertotal) ---------------------------'

		      do i=1,n
			  DX(i) = X(i) - X0(i)
		      enddo


! 		      !geração do modelo físico a partir dos parâmetros de otimização
		      call modelo2d(X)

		      call BOTH(fobjetivo,Gradfobjetivo,X)
		      write(*,*) 'erro relativo (misfit) : ', errorelat*100.0d0,'%'

		      if(fobjetivo.lt.bestFobjetivo)then
			  do i=1,N
			      bestX(i) = X(i)
			  enddo
			  bestFobjetivo = Fobjetivo 
			  bestMisfit =Misfit 
		      endif

!  		      call LBFGS(N,Msave,iter,X,Fobjetivo,Gradfobjetivo,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAGlbfgs)
!  		      call LBFGSmodif(N,Msave,iter,X,Fobjetivo,Gradfobjetivo,DIAGCO,DIAG,IPRINT,EPS,XTOL,&
! 					&  Qlbfgs,RHOlbfgs,ALPHAlbfgs,Slbfgs,Ylbfgs,IFLAGlbfgs)

!  		      call LAMBDALBFGS(N,Msave,iter,X,Fobjetivo,Misfit,Treg,&
! 					&  Gradfobjetivo,gradmisfit,gradTreg,lambdan,&
! 					&  DIAGCO,DIAG,IPRINT,EPS,XTOL,Qlbfgs,Q1lbfgs,Q2lbfgs,RHOlbfgs,  &
! 					&  ALPHAlbfgs,Slbfgs,Ylbfgs,Y1lbfgs,Y2lbfgs,X0,IFLAGlbfgs)

 		      call lmLAMBDALBFGS(N,Msave,iter,X,Fobjetivo,Misfit,Treg,&
					&  Gradfobjetivo,gradmisfit,gradTreg,lambdan,&
					&  DIAGCO,DIAG,IPRINT,EPS,XTOL,Qlbfgs,Q1lbfgs,Q2lbfgs,RHOlbfgs,  &
					&  ALPHAlbfgs,Slbfgs,Ylbfgs,Y1lbfgs,Y2lbfgs,X0,IFLAGlbfgs)
! 





		      write(*,*)'---------------------------------------------------------------------------------------------------------'

		      if(iter.gt.iterant)itertotal=itertotal+1

		      
	      
		      IF(IFLAGlbfgs.LE.0) exit

		      ICALL=ICALL + 1
 
! 		      IF(abs(Fobjetivo-FobjetivOLD).le.10d-2*FobjetivOLD )exit ! para evitar que muitas avaliacoes sejam feitas com a convergencia estagnada
		      IF(errorelat .lt. tolerro)exit
		      IF(ICALL.GT.icallmax) exit ! We allow at most 2000 evaluations of F and G
		      IF(itertotal.GT.maxiterLBFGS) exit ! We allow at most maxiterLBFGS iterations
! 		      if(iter.eq.msave)then
! 			  write(*,*) 'recomecando a atualizacao de H'
! 			  go to 15
! 		      endif
		      FobjetivOLD = Fobjetivo
		      MisfitOLD = Misfit
		  ENDDO !fim das iterações do L-BFGS

! 		  errorelat=dsqrt(bestMisfit)
! 		  write(*,*)'Novo lambdanaux',lambdanaux
! 		  if(kreg.ne.0 .and. lambdanaux.gt.0.2d0 .and. errorelat.gt.tolerro)then
! 		      lambdanaux = 0.8*lambdanaux
! 		      write(*,*)'Novo lambdanaux',lambdanaux
! 		      do i=1,N
! 			  X(i) = bestX(i)
! 		      enddo
! 		      go to 15
! 		  endif
		  


		  deallocate(Qlbfgs)
		  deallocate(Q1lbfgs)
		  deallocate(Q2lbfgs)
		  deallocate(RHOlbfgs)
		  deallocate(ALPHAlbfgs)
		  deallocate(Slbfgs)
		  deallocate(Ylbfgs)
		  deallocate(Y1lbfgs)
		  deallocate(Y2lbfgs)

		  deallocate(DIAG)
		  deallocate(W)
! 		  deallocate(G)
! 		  deallocate(G1)
! 		  deallocate(G2)

	      endif! seleção do método

! 	      call BOTH(fobjetivo,Gradfobjetivo,X)!<só para plotar a solucao final


! ! 	      call modelo2d(X)!<só para plotar a solucao final
! 	      call systempsimage(c0,nx0,nz0,deltaX0,deltaZ0,deltaZ0/2.0d0,deltaX0/2.0d0,'veloX',&
! 		& 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
! 		& wbox,Lz0/4.d0,Lx0/4.d0,&
! 		& 0.0d0,Lz0+(Npml-1)*delta,-(Npml-1)*delta,Lx0+(Npml-1)*delta,1)!modelo do pi, sem a pml
! ! 	      write(*,*)'______ TECLE ENTER PARA CONTINUAR'
! ! 	      read(*,*)
	      call systempsimage(c,nx,nz,delta,delta,0.d0,-(Npml-1)*delta,'veloX',&
		& 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
		& wbox,Lz0/4.d0,Lx0/4.d0, &
		& 0.0d0,Lz0,0.0d0,Lx0,1)!modelo do pd sem PML

	      write(unitX,*)'# kw: ',kw,': ',omegas(kw)/(2.d0*Pi),'Hz',' X   final: '
	      do i=1,nx0*nz0
		  write(unitX,*)X(i)
	      enddo


! 	      call systempsmovieSOL(nx,nz,delta,delta,0.0d0,-dfloat(Npml-1)*delta,0.0d0,1.d0,&
! 		  & filename3,'                    ','profundidade (m)    ',&
! 		  &'distancia (m)       ',wbox,Lz0/4.d0,Lx0/4.d0,0.0d0,Lz0,0.0d0,Lx0,cmin,cmax) !figura c do pd sem PML
! 
! 	      call systempsmovieSOL(nx,nz,delta,delta,0.0d0,-dfloat(Npml-1)*delta,0.0d0,1.d0,&
! 		  & filename3,'                    ','profundidade (m)    ',&
! 		  &'distancia (m)       ',wbox,Lz0/4.d0,Lx0/4.d0,0.0d0,&
! 		  & Lz0+dfloat(Npml-1)*delta,-dfloat(Npml-1)*delta,Lx0+dfloat(Npml-1)*delta,cmin,cmax) !figura c do PD com PML

	      call systempsmovieSOL(nx0,nz0,deltax0,deltaz0,deltaz0/2.0d0,deltax0/2.0d0,0.0d0,1.d0,&
		  & filename3,'                    ','profundidade (m)    ',&
		  &'distancia (m)       ',wbox,Lz0/4.d0,Lx0/4.d0,0.0d0,Lz0,0.0d0,Lx0,cmin,cmax) !figura c0 do pi O(sem PML)


! 	      endif!seleção das freqs

	  enddo ! loop das frequencias


	  if(clust.eq.1)then

	      if(kloop.eq.Nloops)suav=0 !< na última clusterização, não suavizar

	      kw = kw-1

! 	      call clusterizainterat(c0,nx0,nz0)
! 	      if(suav.eq.1)then
! 		  call suaviza(c0,nx0,nz0)
! 	      endif

	      call numclassesIm(Nclasses,X) !< escolhe o numero de classes baseado no erro da nova imagem em relação à original
! 	      call numclassesF(Nclasses,X) !< escolhe o número de classes baseado no menor valor da função objetivo

	      if(Nclasses.gt.0)then !< se Nclasses for igual a zero, não se faz clusterização e X não muda
		  call modelo2d(X)

		  call clusteriza(c0,nx0,nz0,Nclasses)


		  if(suav.eq.1)then
		      call suaviza(c0,nx0,nz0)
		  endif

		  do j=1,nz0
		      do i=1,nx0
			  X((j-1)*nx0+i) = mapCemX(c0(i,j),cmin,cmax) !< altera X de acordo com c0 gerado pela custerização
		      enddo
		  enddo
	      endif

	      call BOTH(fobjetivo,Gradfobjetivo,X)!<só para plotar a solucao final

	  endif !if clust=1


      enddo ! loop do numero de loops




      close(unitX)


      deallocate(c0)

      deallocate(X0)
      deallocate(DX)
      deallocate(X)
      deallocate(bestX)
      deallocate(CX)
      deallocate(xref)
      deallocate(CXR)
      deallocate(Xn)
      deallocate(CXn)
      deallocate(resid)
! 	      deallocate(Jacob)
! 	      deallocate(Hprecond)
      deallocate(deltaR) ! matriz: cada coluna é vetor delta correspondente a um receptor      
      deallocate(JdeX)     ! du/dX
      deallocate(residexpandido)
      deallocate(Gradmisfit) 
      deallocate(GradTreg)
      deallocate(depX)
      deallocate(PreceivBARRA) 
      deallocate(Hvalue)
      deallocate(Hrowind)
      deallocate(Hcolind)
      deallocate(f) !lado direito do sistema do problema direto
      deallocate(rowind)!CCS da matriz de impedancia
      deallocate(colptr)!CCS da matriz de impedancia
      deallocate(values)!valores nao nulos da matriz do sistema do problema direto


      close(unitsol)
      close(unitF)  
      close(unitFgrad)




      return
      end subroutine inversao